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Abstract 

We  compare  several  approaches  to  uncer¬ 
tainty  propagation  that  have  been  used  in  the 
literature  to  formulate  the  uncertainty  in  a  dy¬ 
namical  system  governed  by  ordinary  differen¬ 
tial  equations.  Specifically  we  focus  on  the  evo¬ 
lution  of  probability  density  functions  of  the  as¬ 
sociated  stochastic  processes,  and  discuss  their 
applications  in  different  fields. 

1  Introduction 

Uncertainty  is  ubiquitous  in  physics,  chem¬ 
istry,  bioscience  and  especially  in  the  social  sci¬ 
ence.  It  is  often  classified  into  two  types.  One 
is  epistemic  (or  reducible)  uncertainty,  which 
can  possibly  be  reduced  by  improved  measure¬ 
ments.  The  other  is  aleatory  (or  irreducible) 
uncertainty,  which  is  the  result  of  intrinsic 
variability/stochasticity  of  the  system.  Uncer¬ 
tainty  propagation  through  a  dynamic  system 
has  enjoyed  considerable  research  interests  dur¬ 
ing  the  past  decade  due  to  the  wide  applica¬ 
tions  of  mathematical  models  in  studying  the 
dynamical  behavior  of  the  system. 

In  this  paper,  we  consider  two  types  of  dif¬ 
ferential  equations  in  the  presence  of  uncer¬ 
tainty.  One  is  stochastic  differential  equations 
(SDEs),  a  classification  reserved  for  differential 
equations  driven  by  white  noise.  The  other  is 
random  differential  equations  (RDEs),  a  clas¬ 
sification  used  for  differential  equations  driven 
by  colored  noise.  Specifically,  we  focus  on  the 
evolutions  of  the  probability  density  functions 
of  stochastic  processes  resulting  from  SDEs  and 
RDEs,  and  their  applications  in  different  fields. 


Unless  otherwise  indicated,  a  capital  letter  is 
used  to  denote  a  random  variable,  a  bold  capi¬ 
tal  letter  is  for  a  random  vector,  and  their  corre¬ 
sponding  small  letters  are  for  their  realizations 
throughout. 

2  Stochastic  Differential  Equations 

White  noise  is  widely  used  in  practice  to  ac¬ 
count  for  the  randomness  in  the  inputs  of  sys¬ 
tem  governed  by  differential  equations,  espe¬ 
cially  in  the  case  where  one  has  little  knowl¬ 
edge  about  the  precise  nature  of  the  noise.  In 
this  section,  we  consider  stochastic  differential 
equations  driven  by  Gaussian  white  noise.  Note 
that  stochastic  differential  equations  can  be  in¬ 
terpreted  in  either  Ito  or  Stratonovich  sense, 
and  the  difference  between  their  solutions  dis¬ 
appears  if  the  drift  term  is  corrected  based  on 
the  so-called  Wong-Zakai  theorem  (e.g.,  see  [11, 
Section  3.4]).  Hence,  in  this  section  we  focus 
our  discussions  on  the  Ito  stochastic  differential 
equations. 

Consider  the  following  Ito  stochastic  differen¬ 
tial  equations 

dX(t)  =  g(t,  X(t))dt  +  a(t,  X(t))dW(t), 
x(0)  =  Xo, 

(2.1) 

where  X  =  (Xlf  X2, . . . ,  Xn)T,  g  = 
(<7i,  <?2,  •  •  ■ ,  3n)T  is  a  nonrandom  n-dimensional 
function  of  t  and  x,  a  is  a  nonrandom  n  x  l 
matrix  function  of  t  and  x,  W  (t)  is  an  l- 
dimensional  standard  Wiener  process  indepen¬ 
dent  of  the  random  initial  vector  X0.  In  addi¬ 
tion,  we  assume  that  g  and  cr  satisfy  conditions 
that  guarantee  the  existence  and  uniqueness  of 
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solutions  to  initial  value  problem  (2.1). 

Ito  stochastic  differential  equations  are  of¬ 
ten  found  useful  in  financial  mathematics  as  it 
only  takes  into  account  the  information  about 
the  past  (for  example,  in  modeling  the  stock 
price,  the  only  information  one  has  is  about 
past  events).  Ito  stochastic  differential  equa¬ 
tions  have  also  been  used  in  describing  mech¬ 
anisms  of  climate  variability.  For  example,  it 
was  used  in  [20]  to  describe  the  interactions 
between  the  atmosphere  temperature  and  the 
ocean’s  surface  temperature. 

The  evolution  of  the  probability  density  func¬ 
tion  of  the  solution  to  (2.1)  is  described  in  the 
following  theorem. 

Theorem  2.1.  Assume  thafX.0  has  probability 
density  function  po,  and  X(t)  satisfies  (2.1). 
Then  the  probability  density  function  p  o/X(t) 
satisfies 


3  Random  Differential  Equations 

Colored  noise  has  been  observed  in  many  set¬ 
tings  such  as  astronomy,  acoustic  equipment 
and  electronic  devices.  Hence,  random  differ¬ 
ential  equations  have  gained  considerable  re¬ 
search  interest  in  the  past  decade,  especially 
the  efforts  on  computational  methods  (inter¬ 
ested  readers  can  refer  to  [23]  for  discussions 
on  widely  used  methods  such  as  Monte  Carlo 
methods,  stochastic  Galerkin  methods  and 
probabilistic  (or  stochastic)  collocation  meth¬ 
ods).  In  practice,  one  often  decomposes  and 
truncates  the  colored  noise  by  a  finite  linear 
combination  of  uncorrelated  random  variables 
(e.g.,  using  a  Karhunen-Loeve  expansion)  for 
uncertainty  quantification.  Hence,  in  this  sec¬ 
tion  we  focus  our  discussions  on  those  random 
differential  equations  with  random  inputs  char¬ 
acterized  by  finite  dimensional  random  vectors. 


r \  TL  r\ 

—p(t,  x)  +  far(9k(t,  x)p(t,  x)) 

k=l 


X 


a2 


2  ^  dxkdx 
k,j= i 


(cr(t,  x)<t7  (t,  x))fc_.  p{t,  x) 


j  " 


(2,2) 


with  initial  condition  p( 0,x)  =  p0(x).  Here 
(er(t,  x)<r7  (t,  x))fcj.  denotes  the  ( k,j)th  element 
of  matrix  cr(t,  x)<r7  ( t ,  x). 


Equation  (2.2)  is  often  referred  to  as  Fokker- 
Planck  equation  or  forward  Kolmogorov  equa¬ 
tion,  and  is  important  to  the  fields  of  chemistry 
and  thermodynamics.  Fokker-Planck  equations 
have  also  been  used  in  the  literature  (e.g.,  [3,  8]) 
to  describe  the  population  density  in  a  size- 
structured  population  where  the  growth  pro¬ 
cess  is  a  diffusion  process  satisfying  the  Ito 
stochastic  differential  equation.  There  are  sev¬ 
eral  approaches  that  can  be  used  to  derive  the 
Fokker-Planck  equation.  One  approach  is  to 
employ  expansion  methods  such  as  those  in 
[14,  Section  8.1]  employed  by  Moyal  in  1949. 
Another  method  is  based  on  the  principle  of 
preservation  of  probability  density  functions 
(e.g.,  see  [10]). 


3.1  Differential  Equations  with  Ran¬ 
dom  Initial  Conditions 

In  this  section,  we  consider  the  evolution  of 
probability  density  function  of  the  solution  to 
a  system  of  ordinary  differential  equation  with 
random  initial  conditions 

x  =  g(t,x),  x(0)  =  X0,  (3.1) 

where  x  =  (aq,  X2,  ■  ■  ■ ,  xn)T ,  and  g  = 
(<7i,  (?2,  •  •  • ,  gn)T  is  an  n-dimensional  nonran¬ 
dom  function  of  t  and  x.  Here  X0  is  an  n- 
dimensional  random  vector.  Equation  (3.1) 
is  referred  to  as  a  crypto- deterministic  system 
(e.g.,  see  [11,  14,  18]),  and  has  been  proven 
useful  in  wide  applications  such  as  classical 
statistical  mechanics,  statistical  thermodynam¬ 
ics,  kinetic  theory  and  biosciences.  Depending 
on  the  applications,  the  uncertainty  in  initial 
conditions  could  be  classified  as  either  epis- 
ternic  uncertainty  (e.g.,  the  uncertainty  is  due 
to  the  measurement  error)  or  aleatory  uncer¬ 
tainty  (e.g.,  individuals  have  different  initial 
size). 

Theorem  3.1.  [18,  Theorem  6.2.2]  Assume 
that  (3.1)  has  a  mean  square  solution  x(t;  X0). 


Then  the  probability  density  function  of  the  so¬ 
lution  x(t;  X0)  satisfies 

d  n  c) 

—p(t,  x)  +  gfr(9k(t,  x)p(t ,  x))  =  o, 

fc=i  ^ 

(3.2) 

with  initial  condition  p( 0,  x)  =  p0(x),  where  po 
is  the  probability  density  function  oflK0. 

The  above  result  is  often  referred  to  as  Li¬ 
ouville’s  “ theorem ”,  and  the  resulting  equation 
(3.2)  is  often  called  Liouville’s  equation.  Liou- 
ville’s  “theorem”  is  a  key  result  in  statistical 
mechanics  and  statistical  thermodynamics.  It 
is  worth  noting  that  Liouville’s  equation  with 
n  =  3  is  the  well-know  equation  of  continuity  in 
continuum  mechanics.  There  are  several  meth¬ 
ods  that  can  be  used  to  derive  this  equation. 
One  is  from  a  probabilistic  point  of  view  us¬ 
ing  the  concept  of  characteristic  functions  (e.g., 
see  [18,  pl47])  following  the  derivation  of  Kozin 
[12]  in  1961.  Another  method  is  based  on  the 
principle  of  preservation  of  probability  density 
function  [16,  pages  363-364],  In  addition,  we 
see  that  (3.1)  is  a  special  case  of  Ito  stochastic 
differential  equation  (2.1)  with  a  =  0.  Hence, 
if  we  assume  that  g  satisfies  the  conditions  for 
existence  and  uniqueness  of  solutions  of  (2.1), 
then  we  can  use  the  same  arguments  as  those 
for  deriving  the  Fokker-Planck  equation  to  de¬ 
rive  Liouville’s  equation. 

3.1.1  Application  of  Liouville’s  Equa¬ 
tion  in  Population  Dynamics 

We  observe  that  the  Liouvillc  equation  (3.2) 
with  (j\  =  1  is  a  special  case  of  a  model  given 
by  Oster  and  Takahashi  in  1974  [15]  1 


(3.3) 


without  mortality  (d  =  0).  This  model 
is  used  to  describe  the  population  density 

1The  connection  between  the  Liouville  equation  and 
(3.3)  was  recognized  by  Oster  and  Takahashi  in  the  case 
of  zero  mortality  and  constant  growth  rates  gi,  i  =  2, . . . ,  n. 


u(t,x. )  in  a  spatially  homogeneous  population 
where  individuals  are  characterized  by  x  = 
(xi,X2,  ■  ■  ■  ,xn)T  with  Xi  denoting  the  chrono¬ 
logical  age  and  X2,  ■■  ■  ,xn  representing  some 
physiological  variables  such  as  mass,  volume, 
chemical  composition,  and  any  other  quantities 
having  an  influence  on  individual’s  growth  and 
mortality  rates,  and  g(t,x )  and  d{t ,  x)  are  the 
growth  and  mortality  rates  of  individuals  with 
chronological  age  X\  and  physiological  variables 
X2,  ■■■  ,xn  at  time  f,  respectively. 

Equation  (3.3)  includes  a  number  of  well- 
known  structured  population  models  as  special 
cases.  Specifically,  the  model  with  only  chrono¬ 
logical  age  involved,  that  is, 

d  8 

— u(t,  x)  +  —u(t,  x)  +  d(t,  x)u(t,  x)  =  0 
ot  ox 

is  the  age-structured  population  model  given  by 
Mckendrick  [13]  in  1926  and  Von  Foerster  [21] 
in  1959.  Equation  (3.3)  with  n  =  2, 

d  d  d 

— u(tx)  +  —u(t,x)  +  —(g2{t,x)u(t,x)) 

+d(t,  x)u(t,  x)  =  0, 

(3.4) 

is  the  age-size  structured  population  model  de¬ 
veloped  by  Sinko  and  Streifer  [17]  in  1967  (with 
X2  being  some  physiological  variable),  and  it  is 
also  the  model  given  by  Bell  and  Anderson  [9] 
in  1967  for  cell  populations  (with  x2  denoting 
the  volume  of  the  cell).  Equation  (3.4)  without 
the  second  term,  that  is, 

^u(t,x)  +  ft(g(t,x)u(t,x))  (35) 
+d(t,  x)u(t,  x)  =  0, 

is  often  referred  to  as  the  Sinko- Streifer  model 
or  classical  size-structured  population  model. 
Here  x  is  the  structure  variable,  which  may 
represent  weight,  length,  volume,  chronological 
age  (equation  (3.5)  becomes  an  age-structured 
population  model  in  this  case),  caloric  content, 
maturity,  etc,  depending  on  the  applications. 

Based  on  the  above  discussions,  we  see  that 
all  the  linear  structured  population  models  pre¬ 
sented  in  this  section  can  be  associated  with 
some  stochastic  process,  which  is  obtained  due 
to  the  variability  in  the  initial  size  of  individu¬ 
als  in  the  population. 


3.2  Differential  Equations  with  Ran¬ 
dom  Model  Parameters  and  Ran¬ 
dom  Initial  Conditions 

In  this  section,  we  consider  the  system  of  ran¬ 
dom  differential  equations 

x  =  g(t,  x;  Z),  x(0)  =  X0,  (3.6) 

where  x  =  (x1;  x2, . . . ,  xn)T,  g  = 

(<?!,  g2, . . . ,  gn)T  is  an  n- dimensional  nonran¬ 
dom  function  of  t  and  x,  Z  is  an  m- dimensional 
random  vector  with  values  in  Z,  and  Xo  is  an 
n-dimensional  random  vector. 

Equation  (3.6)  has  wide  applications  in 
physics  and  bioscience  to  account  for  hetero¬ 
geneity  of  the  material  or  the  variation  between 
individuals  (and  the  uncertainty  is  aleatory  in 
this  case).  For  example,  a  system  of  ordinary 
differential  equations  are  often  used  to  model 
the  HIV  dynamics  for  any  specific  individual, 
but  the  model  parameters  vary  across  the  pa¬ 
tients  as  the  clinical  data  shows  a  great  deal  of 
variability  among  HIV  patients  (e.g.,  see  [2,  7]). 
Equation  (3.6)  can  also  be  used  to  describe  the 
uncertainty  arising  from  the  parameter  estima¬ 
tion  error,  and  in  this  case  the  uncertainty  is 
epistemic. 

3.2.1  Evolution  of  Joint  Proba¬ 
bility  Density  Function  of 

(x(t;  X0,  Z),  Z)T 

In  this  section,  we  consider  the  evo¬ 
lution  of  probability  density  functions  of 
(x(t;  X0,  Z),  Z)t.  Let  x  =  (x,  z)7 .  Then  (3.6) 
can  be  rewritten  as 


x  =  g (t,x),  x(0)  =  (X0,  Z)T .  (3.7) 


Then  the  joint  probability  density  function 
y?x,z  (i,  x,  z)  o/x(i;X0,  Z)  and  Z  satisfies 


d_ 

dt 


V^x.z  (t,  x,  z) 


,x;  z)<px>z(t,x.,z))  =  0, 


(3.8) 


with  initial  condition  </?x,z(0,  x,  z)  = 
cdx.z(x’z)>  where  Px.z  * s  the  joint  probability 
density  function  o/X0  and  Z. 


Observe  that  (3.7)  is  a  special  case  of  stochas¬ 
tic  differential  equations  with  diffusion  coeffi¬ 
cient  being  zeros.  Hence,  we  can  also  use  the 
same  arguments  for  deriving  the  Fokker-Planck 
equation  to  derive  (3.8). 

Given  the  joint  probability  density  function 
(p (t,  x,  z),  we  can  obtain  the  probability  density 
function  p(t,  x)  for  x(t;  X0,  Z)  given  by 


p(t,  x)  =  /  <p(t,  x,  z)dz. 


(3.9) 


3.2.2  Evolution  of  Conditional  Prob¬ 
ability  Density  Function  of 
x(t;  X0,  Z)  Given  the  Value  z  of  Z 

In  this  section,  we  derive  the  conditional 
probability  density  function  of  x(i;  X0,  Z)  given 
the  value  z  of  Z.  Note  that  the  conditional 
probability  density  function  ip{t,  x;  z)  is 


z) 


yx,z(i,x,z) 

Pz(z) 


(3.10) 


where  <^z(z)  denotes  the  probability  density 
function  of  Z.  Hence,  by  (3.8)  and  (3.10)  we 
find  ip(t,  x;  z)  satisfies 


Here  g(t,  x)  = 


where  0m  is  an 


g(b  x; z) 

0m 

m- dimensional  column  vector  with  all  the  el¬ 
ements  being  zeros.  Hence,  Theorem  3.1  can 
be  applied  to  (3.7)  to  obtain  the  following  re¬ 
sult  on  the  evolution  of  joint  probability  density 
functions  of  the  solution  to  (3.7). 


Theorem  3.2.  [19,  Theorem  1]  Assume  that 
(3.6)  has  a  mean  square  solution  x(t;X0,Z). 


d  d 

TT <p(t,  x;  z)  +  Q  —  {9k(t, x;  z )<p(t,  x;  z))  =  0, 

k= 1  °  k 

(3.H) 

with  initial  condition  </?(0,x;z)  =  <^o(x;z), 
where  <po  (x;  z)  is  the  probability  density  func¬ 
tion  of  initial  condition  X0  given  Z  =  z. 

Observe  that  for  any  given  value  z  of  Z  sys¬ 
tem  (3.6)  is  crypto-deterministic.  Hence,  by 
using  Liouville’s  equation  (3.2)  the  probabil¬ 
ity  density  function  ip(t,  x;  z)  for  the  solution 


x(f;  X0,  z)  to  (3.6)  with  given  value  z  of  Z  sat¬ 
isfies  (3.11).  Thus,  we  can  derive  evolution 
of  the  conditional  probability  density  function 
</?(£,  x;  z)  directly  from  Liouville’s  equation. 

Given  the  conditional  probability  density 
function  (p{t,  x;  z),  the  probability  density  func¬ 
tion  p(t,  x)  of  x(t;  X0,  Z)  is  given  by 

p(t,x)=  [  (p(t,  x;  z)<pz(z)dz.  (3.12) 
Ji 

We  remark  that  the  scalar  case  of  (3.11) 
along  with  some  proper  boundary  and  ini¬ 
tial  conditions  along  with  (3.12)  have  been 
used  to  describe  the  population  density  in 
a  size-structured  population,  where  the  en¬ 
tire  population  is  partitioned  into  subpopula¬ 
tions  with  each  subpopulation  having  its  own 
growth  rate  g(t,x;  z).  This  model,  referred 
to  as  growth  rate  distributed  size-structured 
(GRDSS)  population  model,  was  first  formu¬ 
lated  in  [1,  6]  in  1986,  and  has  been  success¬ 
fully  used  to  model  mosquitofish  population  in 
the  rice  fields,  where  the  data  exhibits  both 
bimodality  and  dispersion  in  size  as  time  in¬ 
creases  (e.g.,  see  [1]).  In  addition,  this  model 
was  also  used  to  model  the  early  growth  of 
shrimp  populations,  which  exhibits  a  great  deal 
of  variability  in  size  as  time  evolves  even  though 
the  shrimp  begin  with  approximately  similar 
size  [3,  4],  Based  on  the  above  discussions,  we 
see  that  the  GRDSS  model  can  be  associated 
with  some  stochastic  process,  which  is  obtained 
due  to  the  variability  in  the  individual’s  growth 
rate  and  also  the  variability  in  the  initial  size 
of  individuals  in  the  population.  This  is  made 
precise  in  [5,  8]. 

4  Concluding  Remarks 

In  this  note  we  consider  the  evolution  of  prob¬ 
ability  density  functions  of  the  solutions  of  Ito 
SDEs  and  RDEs  as  well  as  their  applications 
in  different  fields.  Even  though  SDEs  and 
RDEs  are  quite  different  in  nature  (one  driven 
by  white  noise  while  the  other  driven  by  col¬ 
ored  noise),  it  was  shown  in  [8]  that  there  are 
classes  of  Ito  SDEs  with  solutions  having  the 


same  probability  density  functions  at  each  time 
t  as  those  for  the  solutions  of  a  related  corre¬ 
sponding  RDE.  More  precisely,  it  is  shown  how 
to  map  the  SDE  to  a  corresponding  RDE  and 
conversely.  In  addition,  based  on  these  point- 
wise  equivalence  results  it  was  demonstrated 
in  [5]  that  the  solution  for  the  GRDSS  model 
can  be  used  to  approximate  the  solution  of 
Fokker-Planck  equation  (which  is  difficult  to 
solve  when  the  drift  dominates  the  diffusion, 
the  case  of  primary  interest  in  many  cases). 
However,  it  should  be  noted  that  the  stochas¬ 
tic  processes  obtained  by  these  SDEs  and  their 
corresponding  RDEs  are  different  as  their  co- 
variance  functions  are  not  the  same  (see  [8]  for 
more  details). 
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